library(vegan)
meth<-read.csv('100methdescB.csv',header=TRUE, row.names=1)
#transform data (spe abundance is high skewed and large values that are valid but highly influenctial) - doing a log transformation but add 1 since log zero is undefined
meth.log<-log(meth+1)
hist.plots(meth.log)
#select a matrix coefficient - here doing a Bray-Curtis
meth.d<-vegdist(meth.log,"bray")
#perform PCoA. k is the number of primical compoents that should be extracted from this distance matrix (max here is min(col,row)-1), eig return eigenvalues default is false, if add is TRUE a constant is added to each value in dissim matrix so that resulting eigenvalues are non-negative (default here is FALSE)
meth.pcoa<-cmdscale(meth.d,k=35,eig=TRUE,add=TRUE)
#look at it
meth.pcoa
#want to know the % variation explained by any principal component? then divide each eigenvalue by the sum of the eigenvalues across all PCs. Below looking at first 5 PCs
meth.pcoa$eig[1:5]/sum(meth.pcoa$eig)*100
#look for signifiance of each PC by compare eigenvalues to expectations according to broken stick model
plot(meth.pcoa$eig/sum(meth.pcoa$eig)*100-bstick(45)*100,xlab="PC", ylab="Actual-random percent variation explained") 
abline(h=0)
#look at ordination plot of PC and PC2
ordiplot(meth.pcoa, choices=c(1,2), type="text", display="sites", xlab="PC1 (14.3%)", ylab="PC2 (11.6%)")
#or maybe want to look at PC 1 and 3
ordiplot(meth.pcoa, choices=c(1,3), type="text", display="sites", xlab="PC1 (15%)", ylab="PC3 (9%)")
#if interested in looking at 1st 3 PCs simultaneously
pairs(spe.pcoa$points[,1:3],pch=24,col="red")
#Now calculate the PCA loadings (i.e. variable weights)
#here performs a linear regression between each of the original desriptors (species)and the scores from each PC. A permutation test is used to assess statistical significance
vec.meth<-envfit(meth.pcoa$points, meth.log, perm=1000)
#look at values
vec.meth
#look closer at whats in vec.sp
names(vec.meth)
names(vec.meth$vectors)
#now plot the eigenvectors on the ordi plot
ordiplot(meth.pcoa, choices=c(1,2), type="text", display="sites", xlab='PC1 (14.3%)', ylab="PC2 (11.6%)")
#same plot just w/0 labels
ordiplot(meth.pcoa, choices=c(1,2), xlab='PC1 (14.3%)', ylab="PC2 (11.6%)")
plot(vec.meth, p.max=0.01, col="blue")

#just variables
ordiplot(meth.pcoa, choices=c(1,2), xlab='PC1 (18%)', ylab="PC2 (8%)")
plot(vec.meth, p.max=0.01, col="blue")
#perform NMBD on log transformed data distrance is distance coefficient to use (default is bray-curtis), k is the number of ordnation axes to generate, autotransform default is TRUE, will transform the data first (usually square root), trymax is the number of random starts to try, default is 20
meth.nmds<-metaMDS(meth.log,distance='bray', k=2, autotransform=FALSE, trymax=100)
#look at it
meth.nmds
#get a list of objects resulting from the analysis
names(meth.nmds)
#evaluate the number of dimensions with a scree plot
nmds.scree(meth.log, distance='bray', k=10, autotransform=FALSE, trymax=20)
#or try the randomization thing
nmds.monte(meth.log, distance='bray', k=3, autotransform=FALSE, trymax=20)
#how good a job does NMDS do? check by looking at corrlation between the caclulated dissimilarities and the plotted values. Specifically, plot original dissimilarities and euclidean distances in the ordination using stressplot()
stressplot(meth.nmds)
#look at a plot of the NMDS
plot(meth.nmds, type='n')
text(meth.nmds, labels=row.names(meth))
#if we'd like to see how a particular descriptor (in this case, a weird fish abundance) changes with location we can make te symbol size proportional to log abundance
plot(meth.nmds, type='n')
points(meth.nmds,cex=meth.log$X.CG.in.MBD.exon)
#if we want to look at the sample scores, which are the coordinates of the samples in k-dimensional ordination, they are stored in the object names points
meth.nmds$points
#to calculate species loadings (i.e. variable weights) on each derived axis, use envfit with the NMDS scores
vec.meth<-envfit(meth.nmds$points, meth.log, perm=1000)
vec.meth
#NOW these loadings can be plotted
ordiplot(meth.nmds, choices=c(1,2), xlab='Axis 1', ylab='Axis 2',)
plot(vec.meth, p.max=0.01, col='blue')